
// Program: adjustq
// Secondary programs: adjust_simple_q, import_matrix_q
*! version 1.0.3 Isabel Onate 9/24/2018
 
 //
cap program drop adjustq
program define adjustq, rclass byable(onecall)
	version 14.1
	syntax [using] [if], PVAL(name) OUTCome(name) Q(name) ///
						[MATrix(name)]  [EXPort(string)] [SAVEdta(string)] ///
						[BY(name)] [replace] [SHEET(string)] [SHEETREPlace] [SHEETMODify] 
						//[genorder(name)]

	*local by "`_byvars'" // to allow for byable(onecall)
	*di in r "`by'"
	*di in r "`_byvars'"
	// ++++++++++++
	// IMPORT DATA
	// ++++++++++++
	
	// Make sure either data or matrix for input
	if ("`matrix'" != "" & "`using'" != "") {
		di as err "You specified to input data both from a matrix and a datset - these options are muttually exclusive"
		exit 101
	}
	
	// From matrix - only saves rows that are p val and outcome?
	if ("`matrix'" != "") {
		clear
		import_matrix_q, pval(`pval') outcome(`outcome') matrix(`matrix') by(`by')
	}
	
	// From dataset
	if ("`using'" != "") {
		use `using', clear
	}
	
	// +++++++++++++++++++++++++
	// RESTRICT DATA - IF OPTION
	// +++++++++++++++++++++++++
	
	if `"`if'"' != "" {
		qui keep `if'
	}
	
	order `outcome' `pval' `_byvars'
	
	// +++++++++++++++++++++++++++++++
	// ADJUSTMENT ACCORDING TO METHOD
	// +++++++++++++++++++++++++++++++
	// This section uses secondary programs defined bellow
	
	if "`by'" != "" {
		qui levelsof `by', loc(by_groups)
		forvalues number_group = 1/`:word count `by_groups'' {
			loc name_group = "`:word `number_group' of `by_groups''"
			tempfile data_bygroup_`number_group'
			preserve
				// restrict sample - by group string
				cap confirm string var `by'
				if _rc == 0 {
					qui keep if `by' == "`name_group'"
				}
				// restrict sample - by group numeric
				cap confirm numeric var `by'
				if _rc == 0 {
					qui keep if `by' == `name_group'
				}
				adjust_simple_q, pval(`pval') outcome(`outcome') q(`q') // simple method	
				qui save `data_bygroup_`number_group''
			restore	
		}
		if (`: word count `by_groups'' > 1 & !mi(`: word count `by_groups'')) {
			use "`data_bygroup_1'", clear
			forvalues number_group = 2/`:word count `by_groups'' {
				loc name_group = "`:word `number_group' of `by_groups''"
				append using "`data_bygroup_`number_group''"
			}
		}
	}
	
	if "`by'" == "" {
		// simple
		adjust_simple_q, pval(`pval') outcome(`outcome') q(`q') // simple method
		// sharp
		*adjust_sharp_q
	}
	
	
	// +++++++++++
	// EXPORTING
	// +++++++++++

	if `"`export'"' != `""' {
		export excel "`export'", firstrow(var) `replace' sheet("`sheet'") `sheetmodify' `sheetreplace'
	}
	if `"`savedta'"' != `""' {
		save "`savedta'", replace
	}
	
	// ++++++++++++++
	// SRORED RESULTS
	// ++++++++++++++
	
	// matrix for stored results
	qui levelsof `outcome', loc(outcomes_list)
	matrix Q = J(_N, 2, .)
	matrix colnames Q = pvalue q
	matrix rownames Q = `outcomes_list'
	
	loc N = _N
	forvalues row = 1/`N' {
		loc pval_fill = `pval' in `row'
		loc q_fill = `q' in `row'
		matrix Q[`row',1] = `pval_fill'
		matrix Q[`row',2] = `q_fill' 
	}

	return loc N = _N
	return loc outcomes `outcomes_list'
	return matrix Q Q
	
end			



// +++++++++++++++++++++++++++++++
//       SECCONDARY PROGRAMS
// +++++++++++++++++++++++++++++++

// Importing matrices
cap program drop import_matrix_q
program define import_matrix_q 
	version 14.1
	syntax, PVAL(name) OUTcome(name) MATrix(name) [by(name)]
	
	// Error messages
	confirm matrix `matrix' // Error will display as is
	local rownames : rownames `matrix'
	local colnames : colnames `matrix'
	
	if `:list posof "`pval'" in colnames' == 0 {
		di as err "Name for p value is not found in matrix columns"
		exit 101
	}

	clear
	qui gen `pval' = .
	qui gen `outcome' = ""
	qui set obs `: word count `rownames''
	
	local rownames : rownames `matrix'
	forvalues n = 1/`: word count `rownames'' {
		qui replace `outcome' = "`: word `n' of `rownames''" in `n'
		qui replace `pval' = `matrix'[`n',`:list posof "`pval'" in colnames'] in `n'
	}
		
	if "`by'" != "" {
		// Error message
		if `:list posof "`by'" in colnames' == 0 {
			di as err "Name for variable with by groups is not found in matrix columns"
			exit 101
		}
		qui gen `by' = .
		forvalues n = 1/`: word count `rownames'' {
			qui replace `by' = `matrix'[`n',`:list posof "`by'" in colnames'] in `n'
		}
	}
		
end		


// Program for adjustment
cap program drop adjust_simple_q
program define adjust_simple_q
	version 14.1
	syntax, PVAL(varname numeric) OUTcome(varname string) Q(name) 
	//[genorder(name)]
	
	// +++++++++++++++++++++
	// ERROR MESSAGES DATA
	// +++++++++++++++++++++
	
	// Confirm variables with p values and outcomes are different
	if "`pval'" == "`outcome'" {
		di as err "Names for variable with p values and outcomes cant be the same"
		exit 101
	}
	
	// Confirm variable with p value exists and its numeric
	cap confirm numeric var `pval'
	if _rc != 0 {
		di as err "Numeric variable with p values is not found in data"
		exit 101
	}
	
	// Confirm variable with outcomes exists and its a string
	cap confirm string var `outcome'
	if _rc != 0 {
		di as err "String variable with outcome names is not found in data"
		exit 101
	}

	// Make sure p values are not missing
	cap assert !mi(`pval') 
	if _rc != 0 {
		di as err "Variable with p values contains missing values"
		exit 101
	}
	
	// Make sure p values are in the [0,1] interval
	cap assert (`pval' >= 0 & `pval' <= 1) if !mi(`pval')
	if _rc != 0 {
		di as err "Variable with p values contains values outside the [0,1] interval"
		exit 101
	}
	
	// Make sure all outcomes are unique (no duplicates)
	cap isid `outcome'  
	if _rc != 0 {
		di as err "Outcomes are not unique"
		exit 101
	}
	
	// Make sure by variable has non missing values
	cap assert mi(`by') 
	if _rc == 0 {
		di as err "Variable with by group contains only missing values"
		exit 101
	}
	
	// Error for existing temp variables
	foreach var in initial_order fdr_temp reject_temp reject_rank total_rejected rank  {
		cap confirm var ``var'' 
		if _rc != 7 & _rc == 111 {
			di as err "Tempvar `var' has been defined before - used in prgram"
			err 101
		}
		if _rc != 111 & _rc== 0 {
			di as err "Tempvar `var' has been generated before - used in prgram"
			err 101
		}
	}
	// Error for existing variables
	foreach var in `q'_string  {
		cap confirm var `var' 
		if !_rc  {
			di as err "The variable `q' already exists in dataset"
			exit 101
		}
	}
	
	
	// +++++++++++
	// ADJUSTMENT
	// +++++++++++
	
	// temp variables
	tempvar fdr_temp
	tempvar reject_temp
	tempvar reject_rank
	tempvar total_rejected 
	tempvar rank
	tempvar initial_order
	
	
	qui {
		// Value needed for calculations
		loc N = _N
		qui gen `initial_order' = _n 		// initial order for p values to sort again at the end
		qui sort `pval'						// sort p values from smallest to largest
		qui gen `rank' = _n if `pval'!=.	// rank of p values
		qui gen `q' = .						// empty variable for q values
		qui sum `pval'					
		loc count = `r(N)'					// To only inlcude p values that are not missing
					
		// lOOP OVER Q VALUES (find the smallest value at which the null can be rejected)		
		forvalues qval = .999(-.001).000 {
			* Generate qr/N
			qui gen `fdr_temp' = `qval'*`rank'/`count'
			* Generate binary variable checking condition p(r) <= qr/N
			qui gen `reject_temp' = (`fdr_temp' >= `pval') if !mi(`fdr_temp')
			* Generate variable containing p-value ranks for all p-values that meet above condition
			qui gen `reject_rank' = `reject_temp'*`rank'
			* Record the rank of the largest p-value that meets above condition
			qui egen `total_rejected' = max(`reject_rank')
			* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
			qui replace `q' = `qval' if `rank' <= `total_rejected' & !mi(`rank')
			* Remove variables
			qui drop `fdr_temp' `reject_temp' `reject_rank' `total_rejected' 
		}			

		qui replace `q' = 0.999 if mi(`q') & !mi(`pval') // we want the q value to be large but not missing it the condition in the previous loop is not met but the p value is not missing. 			
		sort `initial_order'
		*qui drop `rank'
		*order `outcome' `initial_order' `pval' `q'
		if "`genorder'" != "" {
			gen `genorder' = `initial_order'
		}
	}
	


end	
			
/*			
program define adjust_sharp_q
		
end			
*/			



	
	
	
	
	
	
	
	
	
	
